* Process output of TWFE regressions

	foreach occ in physicians lawyers {
		if "`occ'" == "physicians" {
			use "${output}/twfe/geo_twfe/physicians/IndividFixedEffects_ALL_physicians.dta", clear
		}
		else if "`occ'" == "lawyers" {
			use "${output}/twfe/geo_twfe/lawyers/IndividFixedEffects_lawyers.dta", clear
		}
	
	local N_UR = `=N_UR[1]'
	
	egen aux=mean(logptotinc), by(id)
	replace logptotinc=aux
	keep id cz90 __hdfePLACE__ __hdfeINDIVID__ logptotinc
	duplicates drop 	
	gegen mean_fe_id_`occ' = mean(__hdfeINDIVID__), by(cz90)
	gegen count_fe_id_`occ' = count(__hdfeINDIVID__), by(cz90)	
	gegen logptinc_twfesmpl_`occ' = mean(logptotinc), by(cz90)
	rename __hdfePLACE__ fe_cz_`occ'
	keep cz90 *_fe_id* fe_cz* logptinc_twfesmpl_*
	duplicates drop
	isid cz90 
	gen N_UR_`occ' = `N_UR'
	save "${temp}/cz90_twfe_`occ'.dta", replace
	}
	
	* Add CZ-level characteristics 

	use  "${intermediate_data}/cz_characteristics/CZ1990Characteristics.dta", clear 	
	merge 1:1 cz90 using  "${temp}/cz90_twfe_physicians.dta", assert(master match) nogen
	merge 1:1 cz90 using  "${temp}/cz90_twfe_lawyers.dta", assert(master match) nogen
	save "${intermediate_data}/twfe/geographic_variation_cz90level.dta", replace

* Correlation of CZ fixed effects and CZ characteristics 
	
	// Loop over two samples 
	foreach sample in physicians lawyers {
		frames reset
		
		* Create local list with variables 
		# d ;
		loc vars 
			ruralurban logpopulation pop_over65 pop_density cs_educ_ba frac_middleclass
			
			le_agg_q1 adjtreat age65le count_docs_ally count_docs_2017 pcprate nspecdocs_pc npcpdocs_pc 
			medicaidelig_pc medicareelig_pc puninsured2010 reimb_penroll_adj10 GAF 
			
			mean_fe_id_`sample' med_hhinc2016 cz_price_index_hi a0_laspeyres
			job_growth_1990_2010  median_house_value
			;
		# d cr
		
		
		use "${intermediate_data}/twfe/geographic_variation_cz90level.dta", clear

		*  Prepare variables 
		g greater10physicians = count_docs_2017 > 10 
		qui sum fe_cz_`sample'
		gen treateffect = (fe_cz_`sample' - `r(mean)')
		
		qui sum logptinc_twfesmpl_`sample'
		gen logptinc_demean = logptinc_twfesmpl_`sample' - `r(mean)'
		
		foreach y of varlist `vars' {
			qui su `y'
			replace `y' = (`y'-`r(mean)')/`r(sd)'
		}

		* Regressions 
		eststo clear
		cap log close
		log using "${output}/twfe/geo_twfe/`sample'/correlate_CZFE_CZcharacteristics_`sample'.txt", replace text 
		foreach l of varlist `vars' {
			di "*******************************************************************************"
			di "Regressing demeaned CZ-FEs on `: var label `l'':"
			reg treateffect `l', robust
			qui	eststo trt_`l'
			
			if "`l'" == "med_hhinc2016" {
				di "*******************************************************************************"
				di "Regressing demeaned CZ-FEs on `: var label `l'': (restricting to CZs with > 10 Physicians)"
				reg treateffect `l' if greater10physicians == 1, robust
				qui	eststo trt_`l'pcp10			
			}		
			
			di "*******************************************************************************"
			di "Regressing demeaned log(income) on `: var label `l'':"
			reg logptinc_demean `l', robust
			qui	eststo lhs_`l'
		} 
		log close 
		
		foreach p of varlist `vars' {
			local l`p' : variable label `p'
			if regexm("`l`p''", "(mean)") {
				label var `p' "`=substr(`"`l`p''"', 7, .)'"
			} 
		}

		* Plot coefficients
		# d ; 
		label var ruralurban "Urban->Rural Continuum Code 2013";			
		label var npcpdocs_pc "# of PCPs per 100,000";		
		label var nspecdocs_pc "# of Non-PCPs per 100,000";
		label var medicaidelig_pc "# Medicaid Eligible per 100,000";	
		label var medicareelig_pc "# Medicare Eligible per 100,000";		
		label var pop_over65 "Pop. over Age 65";				
		label var logpopulation "Log Population";
		label var adjtreat "FGW Mortality Trt Effect"; 				
		label var age65le "Age 65 LE";
		label var pcprate "PCP Rate";
		label var med_hhinc2016 "Median HH Income in 2016";
		label var count_docs_ally "Total # of Physicans (All years)";
		label var count_docs_2017 "Total # of Physicans (2017)";
		label var GAF "Baseline Smoothed Medicare GAF";
		label var mean_fe_id_`sample' "Mean Individual FEs";
		label var a0_laspeyres "Laspeyres Price Index";
		label var cz_price_index_hi "Laspeyres Price Index (High Income)";
		
		# d cr 
		
		* Preserve labels
		frame create labels
		frame labels {
			set obs 100
			g label = ""
			g indepvar = ""
			tempfile labels
		} 
		loc i = 1
		foreach l in `vars' {
			local label : variable label `l'
			di "`label'"
			frame labels {
				replace indepvar = "`l'" if _n == `i'
				replace label = "`label'" if _n == `i'
			}
			loc ++i
		}
		
		frame labels {
			drop if indepvar == ""
			save `labels', replace
		} 
		
		* Export estimates
		frame create estimates
		frame change estimates
		set obs 1000
		g depvar = ""
		g indepvar = ""
		g estimate = .
		g stderr = .
		g min95 = .
		g max95 = . 
		g N_UR = .
		
		loc j = 1
		loc i = 1
		foreach l in `vars' {
			
			di "Regressing demeaned CZ-FEs on `l'"
				replace indepvar = "`l'" if _n == `i'
				replace depvar = "Demeaned CZ-FEs" if _n == `i'
				estimates restore trt_`l'
				
				replace estimate = _b[`l'] if _n == `i'
				replace stderr = _se[`l'] if _n == `i'			
				replace min95 = _b[`l'] + invttail(e(df_r), 0.975) * _se[`l'] if _n == `i'
				replace max95 = _b[`l'] - invttail(e(df_r), 0.975) * _se[`l'] if _n == `i'
				replace N_UR = e(N) if _n == `i'
				loc ++i 
				
			di "Regressing demeaned log(income) on `l'"
				replace indepvar = "`l'" if _n == `i'
				replace depvar = "Demeaned Log Earnings" if _n == `i'
							
				estimates restore lhs_`l'
				
				replace estimate = _b[`l'] if _n == `i'
				replace stderr = _se[`l'] if _n == `i'			
				local t_critical = invttail(e(df_r), 0.975)
				replace min95 = _b[`l'] + invttail(e(df_r), 0.975) * _se[`l'] if _n == `i'
				replace max95 = _b[`l'] - invttail(e(df_r), 0.975) * _se[`l'] if _n == `i'
				replace N_UR = e(N) if _n == `i'
				loc ++i 		
				
			loc ++j 
				
		}

		di "Regressing demeaned CZ-FEs on hh income (restricting to CZs with > 10 Physicians)"
		replace indepvar = "med_hhinc2016 (physicians>10)" if _n == `i'
		replace depvar = "Demeaned CZ-FEs" if _n == `i'
		estimates restore trt_med_hhinc2016pcp10				
		replace estimate = _b[med_hhinc2016] if _n == `i'
		replace stderr = _se[med_hhinc2016] if _n == `i'			
		replace min95 = _b[med_hhinc2016] + invttail(e(df_r), 0.975) * _se[med_hhinc2016] if _n == `i'
		replace max95 = _b[med_hhinc2016] - invttail(e(df_r), 0.975) * _se[med_hhinc2016] if _n == `i'
		replace N_UR = e(N) if _n == `i'
		drop if depvar == ""
		g sort = _n 
		
		merge m:1 indepvar using `labels', nogen
		replace label = "Median HH Income in 2016 (CZs with > 10 Physicians)" if indepvar == "med_hhinc2016 (physicians>10)" 
		g sample = "`sample'"
		
		tempfile temp`sample'
		save `temp`sample''
	}	
	
	* Export
	append using  `tempphysicians'
	
	drbest_docinc estimate stderr min95 max95  
	drbcount N_UR, gen(N)
	order N_UR, last
	order depvar indepvar label estimate stderr min95 max95 N N_UR
	export delimited using "${mypath}/intermediate_csv/geotwfe09-cz_correlations.csv", replace dataf delim(tab) 		
